home *** CD-ROM | disk | FTP | other *** search
/ Nautilus 1992 July / Nautilus-3-8 / Nautilus-3-8.bin / Tools & Utilities / Techy Stuff / Source ƒ / sky source ƒ / SAT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-16  |  2.8 KB  |  147 lines
  1. #include "sky.h"
  2.  
  3. sat()
  4. {
  5.     double pturbl, pturbb, pturbr;
  6.     double lograd;
  7.     double dele, enom, vnom, nd, sl;
  8.     double q0, v0, t0, j0, s0;
  9.  
  10.     double capj, capn, eye, comg, omg;
  11.     double sb, su, cu, u, b, up;
  12.     double sd, ca, sa;
  13.  
  14.     object = "Saturn      ";
  15.  
  16.     ecc = .0558900 - .000347*capt;
  17.     incl = 2.49256 - .0044*capt;
  18.     node = 112.78364 + .87306*capt;
  19.     argp = 91.08897 + 1.95917*capt;
  20.     mrad = 9.538843;
  21.     anom = 175.47630 + .03345972*eday - .56527*capt;
  22.     motion = 120.4550/3600.;
  23.  
  24.     q0 = 102.28  + 4.092334429*eday;
  25.     q0 *= radian;
  26.     v0 = 212.536 + 1.602126105*eday;
  27.     v0 *= radian;
  28.     t0 = -1.45  + .985604737*eday;
  29.     t0 *= radian;
  30.     j0 = 225.36 + .083086735*eday;
  31.     j0 *= radian;
  32.     s0 = 175.68 + .033455441*eday;
  33.     s0 *= radian;
  34.  
  35.     anom = anom;
  36.     incl *= radian;
  37.     node *= radian;
  38.     argp *= radian;
  39.     anom = fmod(anom,360.)*radian;
  40.  
  41.     enom = anom + ecc*sin(anom);
  42.     do {
  43.         dele = (anom - enom + ecc * sin(enom)) /
  44.             (1. - ecc*cos(enom));
  45.         enom += dele;
  46.     } while(fabs(dele) > 1.e-8);
  47.     vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
  48.         cos(enom/2.));
  49.     rad = mrad*(1. - ecc*cos(enom));
  50.  
  51.     lambda = vnom + argp;
  52.  
  53.     pturbl = 0.;
  54.  
  55.     lambda += pturbl*radsec;
  56.  
  57.     pturbb = 0.;
  58.  
  59.     pturbr = 0.;
  60.  
  61. /*
  62.  *    reduce to the ecliptic
  63.  */
  64.  
  65.     nd = lambda - node;
  66.     lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
  67.  
  68.     sl = sin(incl)*sin(nd) + pturbb*radsec;
  69.     beta = atan2(sl, sqrt(1.-sl*sl));
  70.  
  71.     lograd = pturbr*2.30258509;
  72.     rad *= 1. + lograd;
  73.  
  74.  
  75.     lambda -= 1185.*radsec;
  76.     beta -= 51.*radsec;
  77.  
  78.     motion *= radian*mrad*mrad/(rad*rad);
  79.     semi = 83.33;
  80.  
  81. /*
  82.  *    here begins the computation of magnitude
  83.  *    first find the geocentric equatorial coordinates of Saturn
  84.  */
  85.  
  86.     sd = rad*(cos(beta)*sin(lambda)*sin(obliq) +
  87.         sin(beta)*cos(obliq));
  88.     sa = rad*(cos(beta)*sin(lambda)*cos(obliq) -
  89.         sin(beta)*sin(obliq));
  90.     ca = rad*cos(beta)*cos(lambda);
  91.     sd += zms;
  92.     sa += yms;
  93.     ca += xms;
  94.     alpha = atan2(sa,ca);
  95.     delta = atan2(sd,sqrt(sa*sa+ca*ca));
  96.  
  97. /*
  98.  *    here are the necessary elements of Saturn's rings
  99.  *    cf. Exp. Supp. p. 363ff.
  100.  */
  101.  
  102.     capj = 6.9056 - 0.4322*capt;
  103.     capn = 126.3615 + 3.9894*capt + 0.2403*capt2;
  104.     eye = 28.0743 - 0.0128*capt;
  105.     comg = 168.1179 + 1.3936*capt;
  106.     omg = 42.9236 - 2.7390*capt - 0.2344*capt2;
  107.  
  108.     capj *= radian;
  109.     capn *= radian;
  110.     eye *= radian;
  111.     comg *= radian;
  112.     omg *= radian;
  113.  
  114. /*
  115.  *    now find saturnicentric ring-plane coords of the earth
  116.  */
  117.  
  118.     sb = sin(capj)*cos(delta)*sin(alpha-capn) -
  119.         cos(capj)*sin(delta);
  120.     su = cos(capj)*cos(delta)*sin(alpha-capn) +
  121.         sin(capj)*sin(delta);
  122.     cu = cos(delta)*cos(alpha-capn);
  123.     u = atan2(su,cu);
  124.     b = atan2(sb,sqrt(su*su+cu*cu));
  125.  
  126. /*
  127.  *    and then the saturnicentric ring-plane coords of the sun
  128.  */
  129.  
  130.     su = sin(eye)*sin(beta) +
  131.         cos(eye)*cos(beta)*sin(lambda-comg);
  132.     cu = cos(beta)*cos(lambda-comg);
  133.     up = atan2(su,cu);
  134.  
  135. /*
  136.  *    at last, the magnitude
  137.  */
  138.  
  139.  
  140.     sb = sin(b);
  141.     mag = -8.68 +2.52*fabs(up+omg-u)-
  142.         2.60*fabs(sb) + 1.25*(sb*sb);
  143.  
  144.     helio();
  145.     geo();
  146. }
  147.